{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 教程\n",
    "\n",
    "本教程使用人类样本作为示例.\n",
    "\n",
    "假定我们获得了几个人源基因序列,但不知道他们在人类参考基因组中的确切位置.\n",
    "\n",
    "我们使用PxBLAT将序列与参考基因进行比对,并找到所在位置.\n",
    "\n",
    "## 理解FASTA格式\n",
    "\n",
    "FASTA格式是用文本表示核酸或肽段序列的常见格式.\n",
    "\n",
    "一个FASTA文件中,每个序列条目由描述行和序列行组成.描述行由>符号开始.示例如下:\n",
    "\n",
    "```FASTA\n",
    ">sequence 1\n",
    "ATGCTAGCTAGCTAGCTAGCTAGCTA\n",
    "GCTAGCTAGCTAGCTAGCTAGCTAGC\n",
    "TAGCTAGCTAGCTAGCTAGCTAGCTA\n",
    ">sequence 2\n",
    "GCAGCTAGCTAGCTA\n",
    "```\n",
    "\n",
    "- `sequenc 1`为描述行\n",
    "- 紧接着是序列行\n",
    "- 序列行可以分为多行,且单行长度没有严格限制\n",
    "- fasta可以包含多个序列,每个序列用描述行分隔\n",
    "\n",
    "fasta文件可用在许多生物信息学任务中,如:\n",
    "\n",
    "- 序列比对:识别序列的相似性和不同区域\n",
    "- 数据库搜索:在数据库中查找特定序列\n",
    "- 系统发生学:分析序列之间的进化联系\n",
    "\n",
    "## 准备示例数据\n",
    "\n",
    "`test_ref.fa`,`tesst_case1.fa`已下载至本教程中的`data`目录下."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ">chr1\n",
      "NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN\n",
      "taaccctaaccctaaccctaaccctaaccctaaccctaaccctaacccta\n",
      "accctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaac\n",
      "cctaacccaaccctaaccctaaccctaaccctaaccctaaccctaacccc\n",
      "taaccctaaccctaaccctaaccctaacctaaccctaaccctaaccctaa\n",
      "ccctaaccctaaccctaaccctaaccctaacccctaaccctaaccctaaa\n",
      "ccctaaaccctaaccctaaccctaaccctaaccctaaccccaaccccaac\n",
      "cccaaccccaaccccaaccccaaccctaacccctaaccctaaccctaacc\n",
      "ctaccctaaccctaaccctaaccctaaccctaaccctaacccctaacccc\n"
     ]
    }
   ],
   "source": [
    "!head data/test_ref.fa"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "301 data/test_ref.fa\n"
     ]
    }
   ],
   "source": [
    "!wc -l data/test_ref.fa"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ">case1\n",
      "TGAGAGGCATCTGGCCCTCCCTGCGCTGTGCCAGCAGCTTGGAGAACCCA\n",
      "CACTCAATGAACGCAGCACTCCACTACCCAGGAAATGCCTTCCTGCCCTC\n",
      "TCCTCATCCCATCCCTGGGCAGGGGACATGCAACTGTCTACAAGGTGCCA\n",
      "A\n"
     ]
    }
   ],
   "source": [
    "!head data/test_case1.fa"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 将FASTA转化为2bit格式\n",
    "\n",
    "为了进行序列比对,需要将FASTA转化为2bit格式.\n",
    "\n",
    "PxBLAT实现了`fa_to_two_bit()`函数来进行该过程;`two_bit_to_fa()`则可以执行逆过程."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[0;31mSignature:\u001b[0m\n",
      "\u001b[0mfa_to_two_bit\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\u001b[0m\n",
      "\u001b[0;34m\u001b[0m    \u001b[0minFiles\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0;34m'list[str]'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
      "\u001b[0;34m\u001b[0m    \u001b[0moutFile\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0;34m'str'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
      "\u001b[0;34m\u001b[0m    \u001b[0;34m*\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
      "\u001b[0;34m\u001b[0m    \u001b[0mnoMask\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0;34m'bool'\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mFalse\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
      "\u001b[0;34m\u001b[0m    \u001b[0mstripVersion\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0;34m'bool'\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mFalse\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
      "\u001b[0;34m\u001b[0m    \u001b[0mignoreDups\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0;34m'bool'\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mFalse\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
      "\u001b[0;34m\u001b[0m    \u001b[0museLong\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0;34m'bool'\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mFalse\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
      "\u001b[0;34m\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mDocstring:\u001b[0m\n",
      "Convert one or more FASTA files to two-bit format.\n",
      "\n",
      "Args:\n",
      "    inFiles (List[str]): A list of paths to the input FASTA files.\n",
      "    outFile (str): The path to the output two-bit file.\n",
      "    noMask (bool, optional): If True, do not mask the output sequence. Defaults to False.\n",
      "    stripVersion (bool, optional): If True, strip the version number from the sequence IDs. Defaults to False.\n",
      "    ignoreDups (bool, optional): If True, ignore duplicate sequences in the input files. Defaults to False.\n",
      "    useLong (bool, optional): If True, use the long format for the two-bit file. Defaults to False.\n",
      "\n",
      "Returns:\n",
      "    None\n",
      "\n",
      "This function converts one or more input FASTA files to two-bit format and saves the result to the specified output file.\n",
      "The function takes several optional arguments that control the behavior of the conversion process, such as whether to mask\n",
      "the output sequence or strip the version number from the sequence IDs. The function returns None.\n",
      "\n",
      "Example:\n",
      "    >>> fa_to_two_bit(['input.fasta'], 'output.2bit')\n",
      "\u001b[0;31mFile:\u001b[0m      ~/ana/envs/cctbx/lib/python3.11/site-packages/pxblat/toolkit/__init__.py\n",
      "\u001b[0;31mType:\u001b[0m      function"
     ]
    }
   ],
   "source": [
    "from pxblat import fa_to_two_bit\n",
    "\n",
    "fa_to_two_bit?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from pxblat import fa_to_two_bit\n",
    "\n",
    "fa_to_two_bit(\n",
    "    [\"data/test_ref.fa\"],\n",
    "    \"data/test_ref.2bit\",\n",
    "    noMask=False,\n",
    "    stripVersion=False,\n",
    "    ignoreDups=False,\n",
    "    useLong=False,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "test_case1.fa  test_ref.2bit  test_ref.fa\n"
     ]
    }
   ],
   "source": [
    "!ls data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "将在指定路径生成新的`test_ref.2bit`文件.\n",
    "\n",
    "PxBLAT也支持命令行方式.\n",
    "\n",
    "## 进行序列查询\n",
    "\n",
    "PxBLAT主要使用两个类进行序列比对:`pxblat.Server`和`pxblat.Client`.比对顺序为:\n",
    "\n",
    "1. 初始化服务端\n",
    "2. 使用客户端,发送序列到服务端进行比对\n",
    "\n",
    "服务端一般处于*preparing,ready,stop*状态中的一种.在客户端尝试发送序列时,服务端需要处于*ready*状态.\n",
    "\n",
    "以下是启动服务端的几种方法:\n",
    "\n",
    "### 在上下文中启动服务端"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pxblat import Client,Server"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def query_context() -> None:\n",
    "    '''执行序列查询,使用上下文管理服务端的启停.'''\n",
    "    host=\"localhost\"\n",
    "    port=65000\n",
    "    seq_dir=\"data\"\n",
    "    two_bit=\"data/test_ref.2bit\"\n",
    "\n",
    "    client=Client(\n",
    "        host=host,\n",
    "        port=port,\n",
    "        seq_dir=seq_dir,\n",
    "        min_score=20,\n",
    "        min_identity=90,\n",
    "    )\n",
    "\n",
    "    with Server(host,port,two_bit,can_stop=True,step_size=5) as server:\n",
    "        # work() #假设work是你自己进行准备工作的函数\n",
    "        server.wait_ready()\n",
    "        result1=client.query(\"ATCG\")\n",
    "        result2=client.query(\"AtcG\")\n",
    "        result3=client.query(\"data/test_case1.fa\")\n",
    "        result4=client.query([\"ATCG\",\"ATCG\"])\n",
    "        result5=client.query([\"data/test_case1.fa\"])\n",
    "        result6=client.query([\"cgTA\",\"data/test_case1.fa\"])\n",
    "        print(result3[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Program: blat (v.37x1)\n",
      "  Query: case1 (151)\n",
      "         <unknown description>\n",
      " Target: <unknown target>\n",
      "   Hits: ----  -----  ----------------------------------------------------------\n",
      "            #  # HSP  ID + description\n",
      "         ----  -----  ----------------------------------------------------------\n",
      "            0      1  chr1  <unknown description>\n"
     ]
    }
   ],
   "source": [
    "query_context()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`Client.query()`可以接受各种类型的输入,包括:\n",
    "\n",
    "- fasta文件的路径\n",
    "- 上述路径为元素的列表\n",
    "- 大小写不敏感的核酸/蛋白序列字符串\n",
    "- 上述字符串为元素的列表\n",
    "- 上述文件路径和序列字符串混合组成的列表\n",
    "\n",
    "`Client.query()`返回BioPython的`QueryResult`对象.\n",
    "\n",
    "### 以通用模式启动服务端\n",
    "\n",
    "pxblat.Server也可以手动控制启动和关闭."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pxblat import Client, Server\n",
    "\n",
    "def query_general():\n",
    "    '''执行序列查询,手动管理服务端的启停.'''\n",
    "    host = \"localhost\"\n",
    "    port = 65000\n",
    "    seq_dir = \"data\"\n",
    "    two_bit = \"data/test_ref.2bit\"\n",
    "\n",
    "    client = Client(\n",
    "        host=host,\n",
    "        port=port,\n",
    "        seq_dir=seq_dir,\n",
    "        min_score=20,\n",
    "        min_identity=90,\n",
    "    )\n",
    "\n",
    "    server = Server(host, port, two_bit, can_stop=True, step_size=5)\n",
    "    server.start()\n",
    "\n",
    "    # work() 假设是自定义的某些操作函数\n",
    "    server.wait_ready()\n",
    "\n",
    "    result1 = client.query([\"actg\", \"data/test_case1.fa\"])\n",
    "    # another_work() 假设是自定义的某些操作函数\n",
    "    result2 = client.query(\"data/test_case1.fa\")\n",
    "    server.stop()\n",
    "\n",
    "    print(f\"{result1=}\")\n",
    "    print(f\"{result2=}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "result1=[None, QueryResult(id='case1', 1 hits), None]\n",
      "result2=[QueryResult(id='case1', 1 hits), None]\n"
     ]
    }
   ],
   "source": [
    "query_general()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "输出结果中,`None`表示查询序列没有符合要求的比对结果;`QueryReuslt`则提供了比对的详细信息.\n",
    "\n",
    "除了`Server`和`Client`,PxBLAT还实现了`ClientThread`类,可以指定一个线程来处理查询.\n",
    "\n",
    "## 理解查询结果\n",
    "\n",
    "继续以上下文模式为例,稍作调整,以熟悉查询结果."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pxblat import Client, Server\n",
    "\n",
    "def query_context() -> None:\n",
    "    '''执行序列查询,返回查询结果'''\n",
    "    host = \"localhost\"\n",
    "    port = 65000\n",
    "    seq_dir = \"data\"\n",
    "    two_bit = \"data/test_ref.2bit\"\n",
    "    client = Client(\n",
    "    host=host,\n",
    "    port=port,\n",
    "    seq_dir=seq_dir,\n",
    "    min_score=20,\n",
    "    min_identity=90,\n",
    "    )\n",
    "    with Server(host, port, two_bit, can_stop=True, step_size=5) as server:\n",
    "        # work()\n",
    "        server.wait_ready()\n",
    "        result1 = client.query(\"ATCG\")\n",
    "        result2 = client.query(\"AtcG\")\n",
    "        result3 = client.query(\"data/test_case1.fa\")\n",
    "        result4 = client.query([\"ATCG\", \"ATCG\"])\n",
    "        result5 = client.query([\"data/test_case1.fa\"])\n",
    "        result6 = client.query([\"cgTA\", \"data/test_case1.fa\"])\n",
    "    return [result1, result2, result3, result4, result5, result6]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[None],\n",
       " [None],\n",
       " [QueryResult(id='case1', 1 hits), None],\n",
       " [None, None],\n",
       " [QueryResult(id='case1', 1 hits), None],\n",
       " [None, QueryResult(id='case1', 1 hits), None]]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results=query_context()\n",
    "results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[QueryResult(id='case1', 1 hits), None]"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result3=results[2]\n",
    "result3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "QueryResult(id='case1', 1 hits)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result=result3[0]\n",
    "result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Program: blat (v.37x1)\n",
      "  Query: case1 (151)\n",
      "         <unknown description>\n",
      " Target: <unknown target>\n",
      "   Hits: ----  -----  ----------------------------------------------------------\n",
      "            #  # HSP  ID + description\n",
      "         ----  -----  ----------------------------------------------------------\n",
      "            0      1  chr1  <unknown description>\n"
     ]
    }
   ],
   "source": [
    "print(result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[HSP(hit_id='chr1', query_id='case1', 1 fragments)]"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result.hsps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "HSP(hit_id='chr1', query_id='case1', 1 fragments)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "top_hsp=result.hsps[0]\n",
    "top_hsp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "      Query: case1 <unknown description>\n",
      "        Hit: chr1 <unknown description>\n",
      "Query range: [0:151] (1)\n",
      "  Hit range: [12699:12850] (1)\n",
      "Quick stats: evalue ?; bitscore ?\n",
      "  Fragments: 1 (? columns)\n"
     ]
    }
   ],
   "source": [
    "print(top_hsp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "`top_hsp.query_id`:case1,指的是查询成功的序列的标识\n",
      "\n",
      "`top_hsp.query_range`:[(0, 151)),指的是查询成功的序列的起始位置和结束位置,为左闭右开区间.\n",
      "\n",
      "`top_hsp.query_span`:151,指的是查询成功的序列的片段长度.\n",
      "\n",
      "`top_hsp.query_start`:0,指的是查询成功的序列的开始位置.\n",
      "\n",
      "`top_hsp.query_end`:151,指的是查询成功的序列的结束位置.\n",
      "\n",
      "`top_hsp.query_strand`:1,指的是查询成功的序列的查询方向.\n",
      "\n",
      "以`hit`开始的属性与上述属性类似,不过返回的是参考序列的信息.\n",
      "如`top_hsp.hit_range`:[(12699, 12850)),指的是查询对应参考序列的起始位置和结束位置.\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(f'''\n",
    "`top_hsp.query_id`:{top_hsp.query_id},指的是查询成功的序列的标识\\n\n",
    "`top_hsp.query_range`:[{top_hsp.query_range}),指的是查询成功的序列的起始位置和结束位置,为左闭右开区间.\\n\n",
    "`top_hsp.query_span`:{top_hsp.query_span},指的是查询成功的序列的片段长度.\\n\n",
    "`top_hsp.query_start`:{top_hsp.query_start},指的是查询成功的序列的开始位置.\\n\n",
    "`top_hsp.query_end`:{top_hsp.query_end},指的是查询成功的序列的结束位置.\\n\n",
    "`top_hsp.query_strand`:{top_hsp.query_strand},指的是查询成功的序列的查询方向.\\n\n",
    "以`hit`开始的属性与上述属性类似,不过返回的是参考序列的信息.\n",
    "如`top_hsp.hit_range`:[{top_hsp.hit_range}),指的是查询对应参考序列的起始位置和结束位置.\\n\n",
    "'''\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "top_hsp对象的所有公开属性如下:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['aln',\n",
       " 'aln_all',\n",
       " 'aln_annotation',\n",
       " 'aln_annotation_all',\n",
       " 'aln_span',\n",
       " 'fragment',\n",
       " 'fragments',\n",
       " 'gap_num',\n",
       " 'gapopen_num',\n",
       " 'hit',\n",
       " 'hit_all',\n",
       " 'hit_description',\n",
       " 'hit_end',\n",
       " 'hit_end_all',\n",
       " 'hit_features',\n",
       " 'hit_features_all',\n",
       " 'hit_frame',\n",
       " 'hit_frame_all',\n",
       " 'hit_gap_num',\n",
       " 'hit_gapopen_num',\n",
       " 'hit_id',\n",
       " 'hit_inter_ranges',\n",
       " 'hit_inter_spans',\n",
       " 'hit_range',\n",
       " 'hit_range_all',\n",
       " 'hit_span',\n",
       " 'hit_span_all',\n",
       " 'hit_start',\n",
       " 'hit_start_all',\n",
       " 'hit_strand',\n",
       " 'hit_strand_all',\n",
       " 'ident_num',\n",
       " 'ident_pct',\n",
       " 'is_fragmented',\n",
       " 'match_num',\n",
       " 'match_rep_num',\n",
       " 'mismatch_num',\n",
       " 'molecule_type',\n",
       " 'n_num',\n",
       " 'output_index',\n",
       " 'query',\n",
       " 'query_all',\n",
       " 'query_description',\n",
       " 'query_end',\n",
       " 'query_end_all',\n",
       " 'query_features',\n",
       " 'query_features_all',\n",
       " 'query_frame',\n",
       " 'query_frame_all',\n",
       " 'query_gap_num',\n",
       " 'query_gapopen_num',\n",
       " 'query_id',\n",
       " 'query_inter_ranges',\n",
       " 'query_inter_spans',\n",
       " 'query_is_protein',\n",
       " 'query_range',\n",
       " 'query_range_all',\n",
       " 'query_span',\n",
       " 'query_span_all',\n",
       " 'query_start',\n",
       " 'query_start_all',\n",
       " 'query_strand',\n",
       " 'query_strand_all',\n",
       " 'score']"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[att for att in dir(top_hsp) if not att.startswith('_')]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## PxBLAT 与 BLAT 的 API 比较\n",
    "\n",
    "PxBLAT 对照 BLAT,定义了一些 API,重要的 API 包括:\n",
    "\n",
    "| PxBLAT          | BLAT       |\n",
    "| --------------- | ---------- |\n",
    "| Client          | gfClient   |\n",
    "| Server          | gfServer   |\n",
    "| two_bit_to_fa() | towBitToFa |\n",
    "| fa_to_two_bit() | faToTwoBit |\n",
    "\n",
    "## 实现自己的blat代码\n",
    "\n",
    "以下是对示例的进一步修改,从而实现通用的查询函数.你可以自行修改,让它更符合你的需求."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Program: blat (v.37x1)\n",
      "  Query: TGAGA_151 (151)\n",
      "         <unknown description>\n",
      " Target: <unknown target>\n",
      "   Hits: ----  -----  ----------------------------------------------------------\n",
      "            #  # HSP  ID + description\n",
      "         ----  -----  ----------------------------------------------------------\n",
      "            0      1  chr1  <unknown description>\n",
      "      Query: TGAGA_151 <unknown description>\n",
      "        Hit: chr1 <unknown description>\n",
      "Query range: [0:151] (1)\n",
      "  Hit range: [12699:12850] (1)\n",
      "Quick stats: evalue ?; bitscore ?\n",
      "  Fragments: 1 (? columns)\n"
     ]
    }
   ],
   "source": [
    "from pathlib import Path\n",
    "from pxblat import Client, Server, fa_to_two_bit\n",
    "\n",
    "def check_path(path):\n",
    "    \"\"\"Check that the path exists.\"\"\"\n",
    "    if not Path(path).exists():\n",
    "        raise FileNotFoundError(path)\n",
    "\n",
    "def blat(seqs, ref, host, port):\n",
    "    \"\"\"Align sequences to a reference genome.\"\"\"\n",
    "    ref = Path(ref)\n",
    "    check_path(ref)\n",
    "    two_bit = ref.with_suffix(\".2bit\")\n",
    "\n",
    "    if not two_bit.exists():\n",
    "        fa_to_two_bit(\n",
    "            [ref.as_posix()], two_bit.as_posix(), noMask=False, stripVersion=False, ignoreDups=False, useLong=False\n",
    "        )\n",
    "\n",
    "    client = Client(\n",
    "        host,\n",
    "        port,\n",
    "        seq_dir=ref.parent.as_posix(),\n",
    "        min_score=20,\n",
    "        min_identity=90,\n",
    "    )\n",
    "\n",
    "    server = Server(host, port, two_bit.as_posix())\n",
    "\n",
    "    with server:\n",
    "        server.wait_ready()\n",
    "        return client.query(seqs)\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    # The sequence to align is `test_case1.fa`\n",
    "    seqs = [\n",
    "        \"TGAGAGGCATCTGGCCCTCCCTGCGCTGTGCCAGCAGCTTGGAGAACCCACACTCAATGAACGCAGCACTCCACTACCCAGGAAATGCCTTCCTGCCCTCTCCTCATCCCATCCCTGGGCAGGGGACATGCAACTGTCTACAAGGTGCCAA\"\n",
    "    ]\n",
    "\n",
    "    results = blat(seqs, \"data/test_ref.fa\", \"localhost\", 65000)\n",
    "    print(results[0])\n",
    "\n",
    "    for hsp in results[0].hsps:\n",
    "        print(hsp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "!rm data/test_ref.2bit"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
